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The  response  of  a truncated  conical  shell  to  a 
dynamic  axial  force  is  studied  by  means  of  a large  deflec- 
tion analysis  which  includes  the  concept  of  initial  imper- 
fections . 

The  derived  governing  nonlinear  partial  differ- 
ential equations  are  reduced  to  nonlinear  ordinary  differ- 
ential equations  of  a lower-order  by  a first-approximation 
and  the  Galerkin  method.  For  cones  of  moderate  length, 
the  solutions  are  found  to  be  dependent  on  one  physical 
and  three  geometric  parameters. 

Numerical  solutions  to  these  equations  were 
obtained  by  the  Runge-Kutta  algorithm  and  studies  conducted 
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to  determine  the  effect  of  varying  the  pertinent  parameters. 
It  was  found  that  the  value  of  the  load  parameter  in  the 
region  of  instability,  and  the  number  of  longitudinal  waves, 
increase  significantly  for  increased  rates  of  loading.  The 
results  of  varying  the  geometric  parameters  are  in  agree- 
ment with  trends  indicated  in  experimental  data,  but  not 
predicted  by  previous  theoretical  analyses. 
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CHAPTER  I 


INTRODUCTION 

Thin-walled  truncated  conical  shells  are  playing 
an  ever-increasing  role  in  aero-  and  hydrospace  struc- 
tures, They  act  as  interstage  units,  aero-  and  hydro- 
dynamic  fairings,  space  capsules,  nozzles,  and  reducer 
sections  in  submarine  hulls.  In  order  to  design  such 
structures,  a thorough  knowledge  of  their  response  to 
applied  loads  is  required. 

Considerable  work  in  the  past  has  been  devoted  to 
the  stability  analysis  of  thin  cylindrical  shells,  another 
important  space  structure.  In  this  work,  it  was  found  that 
the  classical  small-deflection  analysis  was  deficient. 

This  led  to  a nonlinear  theory  based  on  finite  deflections, 
and  to  the  concept  of  initial  imperfections.  Recently, 
attention  has  been  focused  on  the  behavior  of  thin  cylin- 
drical shells  subject  to  dynamic  loadings. 

The  stability  analysis  of  conical  shells  has  been 
exclusively  for  static  loadings,  and  has  yet  to  take  in 
the  concept  of  initial  imperfections.  It  is  this  author's 
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intention  to  analyse  the  dynamical  problem;  and  thus,  to 
solve,  as  a limiting  case,  the  static  stability  problem 
with  initial  imperfections. 

1 . 1 Historical  Background 

An  excellent  review  and  discussion  of  past  investi- 
gations of  the  stability  of  circular  cylindrical  shells 
is  available  in  the  article  by  Thielemann  (1) . Also, 
an  equally  good  review  of  the  past  work  in  conical  shells 
is  available  in  an  article  by  Seide  (2) . In  all  the  above 
works,  the  analysis  has  been  static  in  nature. 

Recently,  because  of  the  increased  use  of  shells 
as  aerospace  structures,  their  stability  under  dynamic 
loads  has  been  investigated.  Vol'mir  (3)  was  the  first  to 
approach  this  problem,  using  a finite-deflection  theory. 
This  investigation  was  limited  to  a cylindrical  panel 
subject  to  a rapidly  applied  compressive  load  along  the 
generators.  Typical  response  curves  were  shown.  Agamirov 
and  Volmir  (4)  extended  this  analysis  to  closed  cylindrical 
shells  under  rapidly  applied  hydrostatic  and  axial  loads. 

In  this  and  subsequent  papers,  it  is  assumed  that  buckling 
occurs  with  the  formation  of  a large  number  of  circumfer- 
ential waves;  thus  allowing  the  nonlinear  theory  of  shallow 
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shells  to  be  employed.  Initial  imperfections  coinciding 
in  form  to  the  assumed  deflection  are  used.  By  applying 
restrictions  on  the  time-dependent  part  of  the  radial 
displacement  of  the  shell,  the  problem  is  reduced  to  a 
one  degree  of  freedom  system.  Typical  response  curves 
for  the  hydrostatic  case  are  shown  and  discussed. 
Kadashevich  and  Pertsev  (5)  improved  upon  the  above  by 
relaxing  the  restriction  on  the  deflection  function,  and 
including  the  circumferential  stress  induced  by  the  axi- 
symmetric  inertia  force.  At  the  same  time,  Herbert  (6) 
also  improved  upon  (4)  by  using  a more  general  deflection 
function.  Typical  curves  for  the  case  of  constant-rate 
of  end  shortening  are  given  and  discussed.  Coppa  and 
Nash  (7)  use  a similar  analysis  in  their  investigation  of 
dynamic  buckling.  This  report  also  includes  some  experi- 
mental results. 

The  dynamic  buckling  of  shallow  spherical  shells 
was  investigated  by  Humphrey  and  Bodner  (8) , Budiansky 
and  Roth  (9) ; and  more  recently  by  Koval  and  Bhuta  (10) , 
Ho  and  Nash  (11) . 

Since  the  start  of  the  present  investigation, 
another  article  on  dynamic  buckling  of  cylindrical  shells 
(12)  has  appeared.  No  new  theoretical  techniques  were 
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employed;  however,  the  case  of  a step  axial  load  is 
treated  rather  thoroughly. 

1 . 2 Statement  of  the  Problem 

It  is  the  aim  of  the  present  analysis  to:  i)  derive 
the  differential  equations  which  govern  the  dynamic  non- 
linear response  of  a thin-walled  truncated  conical  shell, 
ii)  extend  the  methods  used  in  the  above  articles  to  the 
solution  of  the  problem  of  a conical  shell  subject  to  a 
dynamic  axial  load,  and  iii)  present  numerical  results 
which  indicate  the  effects  of  various  parameters  on 
response . 

We  specify  the  position  of  a point  on  the  middle 
surface  of  the  conical  shell  by  the  distance  s measured 
from  the  cone  vertex  along  a generatrix,  and  the  angle  jo 
between  the  axial  plane  passing  through  the  point,  and 
the  axial  plane  of  the  origin  of  the  coordinates  (see 
Figure  1)  . Thus,  a point  on  the  surface  can  be  repre- 
sented parametrically  as 

Xj  = S'Sinjr«cos^ 

Xj  = s • siniT  • sin^ 


s- cost 
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Figure  1.  Conical  Shell  Under  Axial 
Compression 
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We  assiune  that  there  exist  initial  deviations  from 
the  ideal  geometric  form,  and  that  these  deviations  are, 
at  most,  of  the  same  order  of  magnitude  as  the  thickness 
of  the  shell.  We  shall  then  construct  a theory  for  the 
nonlinear  asymmetric  response  of  a thin-walled  truncated 
conical  shell  to  a time-dependent  axial  load.  The  type 
of  nonlinearity  introduced  is  the  usual  "geometric  non- 
linearity. " 


CHAPTER  II 


ANALYTICAL  FORMULATION 


2 . 1 The  Deformation  and  Elasticity  Relations 

We  shall  assume  that  the  initial  imperfections 
and  the  deformation  resulting  are  such  that  they  divide 
the  shell  into  a number  of  shallow  parts.  This  allows 
us  to  use  Donnell,  Mushtari-Vlasov  strain-displacement 
relations  which  for  a conical  shell  with  initial  imper- 
fections become 
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This  assumption  also  permits  us  to  use,  within  the 
same  degree  of  approximation,  linear  change  of  curvature 
and  twist  relations 
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We  use  the  linear  relations  suggested  by  Novozhilov 

(13)  between  the  forces,  moments,  and  the  strains  of  the 
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and 
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Figure  2. 


Force  and  Moment  Resultants 
Acting  on  a Shell  Element 
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The  above  assumption,  and  the  restriction  placed 
on  the  type  of  nonlinearity  considered,  results  in  a 
theory  which  is  referred  to  in  the  literature  as  the 
nonlinear  theory  of  shallow  shells. 

2 . 2 Hamilton's  Principle  and  the  Resulting  Equations 
of  Motion 

We  shall  obtain  the  governing  differential  equations 
from  Hamilton's  principle  for  conservative  systems  which 
states: 

The  motion  of  the  system  from  time  tj  to  time  tg  is 
such  that  the  line  integral 


I = 


L dt 


where  L=T-V,  is  an  extremum  for  the  path  of  motion. 
That  is,  the  variation  of  the  line  integral  I for  fixed 
and  tg  is  zero.  Symbolically, 


§ 


L dt 


0 
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The  kinetic  energy  of  the  shell  can  be  taken  as 


r'2-n'sinif'  ns, 


T 


St  ) 


where  the  longitudinal  and  circumferential  kinetic  energy 
terms  have  been  neglected,  which  is  compatible  with  the 
assumption  that  the  shell  divides  into  a number  of  shallow 
parts.  However,  since  the  load  is  introduced  into  the  shell 
along  the  generators,  neglect  of  the  longitudinal  inertia 
does  restrict  the  loading  rate. 

The  total  potential  energy  V is  composed  of  the 
strain  energy  Uj_  stored  in  the  elastic  body  due  to  the 
actual  deformations,  plus  the  potential  energy  Ug  of  the 
external  loads.  These  quantities  can  be  expressed  in  the 
form 


^2ffsinr  psg 


1 

2 


( Ng€g  + + Si*  Qp 


0 


s, 


+ Mgkg  + M^k^  + 2Hkg^  ) sdsd^i/ 


(2.4) 


and 
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p2~sin'!f' 


U 


[s  ( NgU  + Ng^v  - QgW 


, „ 3W  ^ 1 3W 

+ Mc3  r:—  + 

S ss  Sf  S 30 


)] 


(2.5) 


Replacing  the  strains  by  the  mid-plane  force  resultants; 
the  moment  resultants,  change  of  curvature,  and  twist  by 
their  displacement  relations;  and  performing  the  indicated 
variation,  yield  the  following  set  of  equations 


) - N.  + 


= 0 


(2.6a) 


3N<o  ^13 
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(2.6b) 
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where 
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plus  a set  of  boundary  conditions,  which  along  with  the 
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differential  equations  will  render  the  integral  I an 
extremum  for  the  path  of  motion.  Neither  these  con- 
ditions nor  their  derivation  will  be  presented  here, 
since  they  have  been  obtained  in  a previous  work  (14) . 

2 . 3 Introduction  of  a Stress  Function 

Equations  (2.6a)  and  (2.6b)  are  identical  to  the 
equilibrium  equations  (in  polar  coordinates)  for  two- 
dimensional  problems  in  the  theory  of  elasticity.  There- 
fore, we  can  reduce  the  number  of  unknowns  by  the  intro- 
duction of  the  Airy  stress  function.  Thus, 


N 


s 


i ( ;^  + 1.  ) 

s Bs  s 


(2.7a) 


N, 


(2.7b) 


S 


S 

Bs 


( 1 M ) 


(2.7c) 


thereby  satisfying  equations  (2.6a)  and  (2.6b)  identi- 
cally. To  yield  a possible  force  distribution,  this 
function  must  insure  that  the  condition  of  compatibility 


is  satisfied. 
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2 . 4 Compatibility  Condition 

The  compatibility  condition  for  shells  is  the 
relation  between  the  strains  and  the  normal  component 
of  the  deformation,  which  must  be  satisfied  in  order  to 
insure  the  existence  of  functions  u and  v connected  with 
the  strain  components  by  equations  (2.1).  This  relation 
is  easily  obtained  from  the  strain-displacement  relations 
by  performing  the  operation 


yielding 


-3  6 s _ ^ 


s-ctniT  ( '3^' 


, 1 a r , -sw,  ,2 

> ■ 2 ® ^ ^ ) 


a 

0s  ■ 
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' ' SS* 


(2.8) 
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Substitution  of  the  inverse  of  equations  (2.3) 
and  the  relations  (2.7)  into  the  above,  and  rearranging, 
we  obtain  the  requirement  on  i to  insure  that  the  com- 
patibility condition  is  satisfied. 


ctnH' 


■ss^ 
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s^ 


t 2 
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S "SS  ^S^J 
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ctnzr 


•a  We 
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^ Wa  ' 2 
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_1  a^Wo 

^ s^  as^  a^^ 


. 2 BK  cfwe  _ _1  . avfe 

s^  a^  asa^  ’ 


1.  awe  a^v^' 
s as  as* 


(2.9) 


Thus,  the  problem  of  the  nonlinear  response  of  a 
conical  shell  to  a time-varying  axial  load  has  been  re- 
duced to  the  solution  of  two  4th-order  nonlinear  differ- 
ential equations  (2.6c)  and  (2.9),  subject  to  appropriate 
boundary  and  initial  conditions. 

Equations  (2.6c)  and  (2.9)  reduce  to  those  de- 
rived by  Schnell  (14)  by  setting  the  initial  deflection 
to  zero  and  neglecting  the  inertia  term.  These  equations 
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may  also  be  obtained  from  those  for  shallow  shells  of 
revolution  given  by  Mushtari  and  Galimov  (15) . 

2 . 5 Limiting  Cases 

It  is  easily  shown  that  the  above  equations  reduce 
to  the  following  well-known  results: 

i)  The  nonlinear  equations  of  motion  for  a 
cylindrical  shell  with  initial  imperfections 
(Donnell ' s equations)  can  be  obtained  by 
letting  the  conicity  approach  zero,  the  dis- 
tance s approach  infinity,  and  ~c't'n/  approach 
some  constant  value  R. 

ii)  The  nonlinear  equations  of  motion  for  a 
circular  plate  with  initial  imperfections 
(Marguerre ' s equations)  can  be  obtained  by 
letting  the  conicity  approach  . 

2 . 6 Transformations  and  Substitutions 
Multiplying  the  operator  A A by  s'^ , we  have 


s^A  A = 
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as*  + 2s'^  -5s* 
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S + S -5>s 


+ 2s“= 


- 2s 


•SSS^ 


■z  + 4 
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which  is  easily  recognized  as  an  operator  of  the  Euler 
type  in  the  variable  s.  Therefore,  it  can  be  trans- 
formed into  a linear  operator  with  constant  coefficients 
by  the  substitution 

2 

s = Sj  e 

Let  the  symbols  ' = and  ‘ = 

The  equation  of  motion  becomes 
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IV 
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+ Ng  ( w,"  - w;'  ) + 2S  ( w;'  - w,  ) 


- s^e“^h 
S|  e pon 


(2.10) 


and  the  compatibility  equation  (2.9)  becomes 


^also  referred  to  as  Cauchy,  homogeneous, 
equidimensional 
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The  above  relations,  and  those  for  strain-displacement 
as  a function  of  z and  , suggest  the  substitutions 


s.^e^^F 


w.  = s.  e ctnlfw, 


Wj,  = s , e ctn5f Wq 
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Then, 


N_  = F + 2F  + F 

O 


= f"  + 3f'  + 2F 


S = - ( F + F ) 


(2.12a) 

(2.12b) 

(2.12c) 


The  equation  of  motion  and  the  compatibility  equation 
become 
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+ Ng  ( + ^'  ) + 2S^,  - s,2e“ 


(2.13) 


f''^  + 4f"'  + 4f“  + 2f"  + 4f'  + 4F  + f 

= Eh  ctn^/j^  [ W|'  - ( 1 + (^  + ^,' + W|  ) ( + W|'  )] 

- - ( 1 + ^;+  ^^ ) ( )]j  (2.14) 

The  main  significance  of  this  substitution  is  the  elimi- 
nation of  the  exponential  function  from  the  compatibility 
equation;  a result  which  is  desirable  for  the  solution 
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technique  to  be  employed.  This  is  the  set  of  differential 
equations  for  which  a solution  will  be  sought. 

2 . 7 Some  Further  Relations 

At  this  point,  we  shall  derive  some  relations  which 
will  be  needed  in  the  ensuing  analysis. 

2.7.1  Periodicity  condition  .-  Since  the  shell  is 
closed  with  respect  to  the  <p  coordinate,  we  have  the  re- 
quirement of  periodicity.  The  condition  that  v be  a 
periodic  function  of  gp  , yields  a relation  which  may  be 
used  to  relate  the  free  parameters  in  an  assumed  solution 
form 


o2TT 


dps  = 0 


(2.15) 


which,  for  our  particular  set  of  variables,  becomes 


^2ffsinJ^ 


ctnif  V d^  = 0 


(2.16) 


Now,  the  strain-displacement  relations  are 
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6 


s 


ctnV  ( u + u ) + -^  ctnV  [ { + 


w 


2 

) 


/ m ,2  1 

- ( W„  + Wo  ) j 


= ctmru  - 


ctn^*'  w + 


ctnir'^  + 


1. 

2 


w,  - tVo  ) 


ctnJi'v' + ctnjfu  + ctn^jf  ^ ) 


-.  % ( % )] 


From  (2.17b),  we  obtain 


ctniTV  = £y>  + 


CtnV  w - 


ctnV  ■ ( 


w,  - 


w‘)  - Ctnjru 


and 


ctny^' = + ctn^a'w'-  ctn^2^  • ( w/w,  - wj^,  ) - ctnyu* 


Adding  yields 


ctn^  ( ^ + ^') 


€.^  + Gqp 


ctn^if 


w + 


w ) 


(2.17a) 


(2.17b) 


(2.17c) 
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[( 


W,  - 


% 


) + 2 ( 


W,  W|  - w,  w^ 


>] 


- ctniT  ( u’  + Q ) 

Substituting  into  the  above,  the  expression  for 
ctnV  ( u' + Q ) obtained  from  (2.17a),  and  the  expression 
for  ctnifv' obtained  from  operating  on  (2.17c),  yields 


ctnirv 


— €.p  + sp  ~ ~ ^s^  ctn^2^  ’ L ^ 2 + W| 


+ w,  + 2^  ) ( ^'+  ^,  ) - ( 2 + wj  + + 2w„  ) ( 


+ Wg  ) + w,  - Wg  j + ctn^u 


Using  relations  (2.3a,  b,  c) , and  requiring  sufficient 
smoothness  in  u,  we  obtain 


p2  TTsin^^ 


( Ehctn^JT  [ 


- Ng  - 


2S  + + u(  N, 


0 


- Ng  - 2S  - Ng 


)]  + ■^  [(  2 + w/+  ^ 


+ 2t$,  ) ( w'  + w,  ) - ( 2 + Wo  + w,  + 2^  ) ( w,' 


+ = 


0 


(2.18) 
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a relation  which  F and  w should  satisfy. 

2.7.2  Unit  end  shortening  .-  We  shall  define  the 
mean  unit  end  shortening  as 


£ 


OSo 


Sg-Sj 


su 

as 


ds 


where  the  bar  indicates  a mean  with  respect  to  the  ^ 
coordinate.  In  terms  of  the  unknown  functions  w and  F, 
we  have 


e = 


s.-s, 


s,e 


ctn^y  • [(  w,'  + w,  f - ( + w^  f] 


Eh 


( f'  + 2F  + F ) 


4/(  f''+  3f'  + 2F 


)] 


dZ 


(2.19) 


where 


With  this  relation  (2.19),  a load  versus  end  shortening 


diagram  can  be  obtained. 


CHAPTER  III 


METHOD  OF  SOLUTION 

Because  of  the  complexity  of  the  governing  differ- 
ential equations  (2.13)  and  (2.14),  recourse  to  an 
approximate  solution  technique  seems  called  for. 

3 . 1 Choice  of  Solution  Form 

We  shall  seek  a solution  in  the  form 

W|  = f^  (t)  sinaasin/3^  + f^  (t)  sinVasin^/^  (3.1) 

where  the  f^ (t)  are  unknown  time  functions.  With  respect 
to  the  space  configuration,  the  first  term  corresponds  to 
that  used  in  the  stability  theory  of  infinitesimal  de- 
flections. The  second  term  reflects  the  preferred  inward 
bulging  of  the  shell  when  the  displacements  become  large. 
The  parameters  cx.  and  are  equal  to  and  ^ respec- 
tively, where  m - number  of  half-waves  along  a generatrix, 
and  n = number  of  full-waves  along  a parallel  circle. 

The  increase  in  wave  length  and  amplitude  of  w as  one 
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moves  in  the  positive  direction  along  a generator  is  in 
agreement  with  the  observed  buckling  pattern.  Also,  for 
the  limiting  case  of  a cylindrical  shell,  this  solution 
form  agrees  with  that  used  by  (4) , (6) , and  (7)  in  their 

analyses  of  the  related  problem. 

The  assumed  displacement  function  satisfies  the 
geometrical  boundary  condition  of  w = 0 at  E.  = 0 and 

Z = Z^.  However,  the  support  condition  at  the  bound- 
aries is  neither  pinned  nor  clamped.  This  rather  arbi- 
trary degree  of  end-fixity  should  be  of  little  signifi- 
cance, provided  the  shell  is  of  moderate  length. 

We  assume  that  only  the  initial  imperfections 
which  have  the  same  form  as  the  solution  are  significant; 
an  assumption  which  is  in  accord  with  the  results  of 
Timoshenko  and  Gere  (16)  for  an  axially  loaded  column. 
Thus, 


= fj^i  sino-asirys^^  + f^^  sino<E  sin^f^  (3.2) 

where  are  determined  from  empirical  formulas  which 
"lump"  all  the  imperfections  into  an  equivalent  geometric 
deviation. 


/ 
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3 . 2 The  Stress  Function 

Introducing  (3.1)  and  (3.2)  into  the  compatibility 
equation  (2.14),  we  obtain  the  particular  solution  by  the 
usual  method  of  undetermined  coefficients.  For  the  homo- 
geneous part  of  the  stress  function,  we  assume  that 

= F„  ( E ) ; thus,  F^  satisfies  the  differential  equation 

f"^  + 4f'"  + 4f"  = 0 

This  equation  yields  the  solution 

^ ) ~ ^1  Cg£e  + CjE  + C^e 

We  set  = 0 since  this  term  does  not  produce  in-plane 
forces.  The  third  term  represents  an  extraneous  root, 
since  the  differential  equation  is  of  higher  order  than 
required  for  the  axisymmetric  problem;  and  therefore,  will 
be  discarded.  Thus, 

F = Ehctn^^  ( J,  sin«2  sin/«;^<  + Jg^costxE  sin/3 
+ Jj  sin2«E  + J^coszocE  + J^cos^;^ 

+ sinZKZcos^;^  + CO  sax  ECO  s^^ 

+ sin3xasin/3j4/ + J,  cos3o(Zsin/3^  + 
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+ J^^sin«£  sin:^^  + J,,  coso<€  sin3/s;^ 

+ J^2.sin3«£sin3/3;^  + J^cossiXE  sin3/S;i/ 
+ J,^sin4o^a  + J,5C0S4«£  + J,4^cos4j/3^ 

+ J,^  sin  40^2  cos  + J,g,cos4*{£cos2)^j4^ 

+ J|q  sin2o<£cos^^  + Jg^coszixE  cos^;^ 
+ Jg|  sin4«£cos4/3^  + J^cos4^£cos4^;i^ 


and  the  Cj_  are  chosen  such  that  force  equilibrium  in  the 
mean  is  satisfied  at  the  boundaries.  Explicitly, 


+ Cj  + C^gZe 


-22 


(3.3) 


where 


Ji  ( f j , f 2 / fji  / f^^  , (X  , ^ ) 


@s  = s 


P + 2tts,  sinlJ'NgCOSiT  - 


(g)s  = Sg 


^S2  p2Tfsin)^ 


J s,  J 0 


0 
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Figure  3.  Equilibrium  of  Shell 


3 . 3 The  Galerkin  Integrals 

We  shall  now  determine  the  unknown  time  functions 
f^(t)  by  the  Galerkin  method.  The  finite  sums  (3.1)  and 
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(3.2)  and  the  stress  function  ordinarily  will  not  satisfy 
equation  (2.13),  but  will  result  in  some  residual  (error 
function)  which  we  will  designate  e^  Thus,  the 

problem  is  to  select  the  f j_  (t)  so  as  to  minimize  e^  (2:,^). 
In  the  Galerkin  method,  one  minimizes  e^  (Z,^)  by  im- 
posing on  it  a set  of  orthogonality  conditions  with  respect 
to  the  space-coordinate  functions.  In  our  case. 


''2Tfsina' 


(z,^)  sincxisin^y/dad)^  = 0 (3.4a) 


p2‘rrsin^' 


(z,y^  ) sin^oCZ:  sin^y/dad^  = 0 (3.4b) 


0 Jo 


These  integrals  are  referred  to  in  the  literature  as  the 
Galerkin  or  Bubnov-Galerkin  integrals. 

Carrying  out  the  indicated  substitutions  and  in- 
tegrations yield,  after  considerable  labor,  the  following 
equations  which  relate  the  unknown  deflection  functions 
with  the  time-varying  axial  load 


d^f,  _ 


^3^1  ^2. 


+ ^8"*” 


(3.5a) 
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( b,+  bjf^+  b^fff^ 

+ ’^8^'+  *’9'^  < *^lo  ’’li  V 3 (3.5b) 

where 


a, 

= 

a,  ( oc  , Z„) 

= 

ag  ( o<  , ^ , Zo) 

^3 

= 

a^  ( o4  / f Z^  / K / f^i 

^4 

= 

a^  ( o(  ,/3  ) 

^5 

= 

a^ioi  ,/3  , Z^,  K) 

a® 

= 

a^  ( oc  f ^ ) 

^7 

= 

a -j  ( oc  , /3  / fg, , f^f) 

^6 

= 

a g ( o<.  , /3  / Zg  1 K , f 

= 

a,  (oc  , Zl„,  K) 

^1 

= 

b,  (oc  ,Z^) 

= 

bg  ( 0^  / > 2,) 

^3 

= 

bg  ( 04  , ^ , E^/  K , fg| 

b4 

= 

b^  ( o4  , ^ , 2.^  / K ) 

^5 

= 

b^  ( «?c  , /3  ) 

= 

b^  ( oc  , /3  ) 

^7 

=: 

b-^  ( 04  , ^ , f p I , f 

= 

bg  ( (X  ,/S) 
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^9  ~ ^9  ( °^  / /3  / K / fjji  / f„2.) 

b,^  = b,„  (<X  K) 

b„  = b„  ( o<-  , /3  , 2.,,  K) 

and  we  have  introduced  the  following  dimensionless 
quantities 

t = is2  t 

R^l 

, , P/S|Tr  sin2«' 

q = V3(l-^')  (-  ) 

K = *^12  (1-A"^)  ( ) ctn^K* 

h 

q can  be  physically  interpreted  as  the  dimensionless 
longitudinal  compressive  stress  at  s,  , normalized  such 
that  for  classical  buckling  (cylinder  or  cone)  q = 1.0. 

The  expression  for  unit  end  shortening  (2.19)  be- 
comes 


£ 


d,  f,+ 


d,f,+ 


where  the  d^  are  given  explicitly  in  the  appendix. 
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3 . 4 Determination  of  Initial  Imperfection 

Initial  imperfections  of  a shell  may  be  the  result 
of;  i)  deviations  from  the  assumed  ideal  geometric  shape, 
ii)  initial  stresses,  iii)  deviations  from  isotropy  of 
the  shell  material,  etc.  Following  Donnell  (17) , we  shall 
assume  that  the  end  result  of  all  the  above  imperfections 
can  be  obtained  by  introducing  an  equivalent  initial  geo- 
metric deviation.  Obviously,  the  effect  of  an  initial 
deflection  depends  not  only  on  its  amplitude,  but  also  on 
the  dimensions  of  the  part  of  the  shell  under  consideration; 
in  our  case  R^,  h,  m,  n.  Therefore,  instead  of  using  an 
arbitrary  constant  value  of  the  amplitude  of  initial  de- 
flection, we  shall  employ  the  concept  of  an  "unevenness" 
factor  introduced  by  Donnell  and  Wan  (18) ; however,  we 
shall  use,  with  a slight  variation,  the  relation  suggested 
by  Loo  (19) . Let 


where  Lg  is  the  length  of  the  first  half-wave  along  a 
generator,  L<jj  is  the  half-wave  length  along  the  smallest 
parallel  circle,  and  Uj,  is  a measure  of  the  unevenness. 


a 
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From  the  above  definitions,  we  have 


1 

r s 9 tn  _ 
[ ( s7  ) “ 1 


and  thus. 


ttR^i  cos^r 


n 


= U J ( 


R 


pi- 


h 


[(  “57  - l]  ctnV' 

1 + WH7  [(  -57  >”  - 1] 


Now, 


f 


01 


R 


h 


yielding. 


01 


= ( 


R 


P\ 


h 


1 

m 


[(  It  )"  - 1] 


1 + - 


n 


(3.6) 


TTsinii* 


[(17)™-  1] 


In  the  numerical  work,  it  was  found  that  for 
equal  to  fractional  values  of  f^^,  the  results  were  in- 
sensitive to  this  parameter.  Therefore,  we  shall  set 
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f 


02. 


.Olf. 


CSI 


(3.7) 


3 . 5 Pertinent  Parameters 

Substitution  of  (3.6)  and  (3.7)  into  (3.5)  renders 
the  equations,  functions  of  the  following  parameters 


( "^  ) / ( ) ,^,Uo,q,Z/,m,n 

Assuming,  for  the  present,  that  the  correct  response  mode 
can  be  determined  from  an  inspection  of  the  results,  and 
that  for  the  shells  of  interest  Uo  and  4/  are  essentially 
constant,  the  pertinent  parameters  become 

^ ¥7  ) ' ( ~ir  ) > t , q 


the  first  three  being  geometric  in  nature;  the  fourth, 
physical . 

3 . 6 Numerical  Solution  Technique 

The  approximate  governing  differential  equations 
(3.5)  can  be  easily  uncoupled  in  the  acceleration  terms, 
resulting  in  the  set  of  equations 
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{t  • f,  , f^) 


i = 1,  2 


These  two  nonlinear  differential  equations  of  the  2nd- 
order  can  be  written  as  a set  of  Ist-order  equations  by 
letting 


df, 

d* 

df^ 


yielding  the  system 


df 


i 


- ^i  (*6;  f|  / f 2^ , f 3 , f^.) 

i = 1,  2,  3,  4 


(3.8) 


subject  to  initial  conditions 


f.  = f 

• 01 


^3 


= fot 


^4  ^oz 


Since  the  system  is  nonautonomous , a phase- space  top- 
ological procedure  to  determine  its  behavior  is 
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impracticable;  and,  thus,  a numerical  technique  will  be 
used.  System  (3.8)  is  amenable  to  the  Runge-Kutta  algo- 
rithm; however,  from  a practical  standpoint,  only  if 
high-speed  electronic  computers  are  available.  Thus,  a 
solution  program  was  written  in  "Fortran  language, " and 
the  numerical  work  was  carried  out  on  the  University  of 
Florida  Computing  Center's  IBM  709  (see  appendix  for  the 
complete  Fortran  program) . 


CHAPTER  IV 


APPLICATION  OF  ANALYSIS 

Although  the  equations  (3.5)  are  valid  for  arbi- 
trary q (ifr)  , in  the  application  of  the  analysis  we  shall 
restrict  ourselves  to  the  case  where  q(si)  is  a linearly- 
increasing  function  of  Tfc,  i.e.  constant-rate  of  axial 
compression.  We  shall  investigate  not  only  the  effect  of 
different  rates  of  compression,  but  we  shall  also  investi- 
gate the  effects  of  the  various  geometrical  parameters  of 
the  problem. 

4 . 1 Selection  of  Response  Mode 

A typical  load  parameter  versus  end  shortening 
response  curve  is  shown  in  Figure  4.  We  see  that  upon 
deviating  from  the  linear  response,  a rapid  increase  in 
end  shortening  occurs;  and  upon  reaching  a maximum  end 
shortening,  the  response  consists  of  nonlinear  oscil- 
lations. We  shall,  however,  not  concern  ourselves  with 
these  oscillations,  since  it  is  highly  improbable  that 
the  shell  is  still  in  an  elastic  state.  Of  greatest 
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Figure  4.  Typical  Response  Curve 
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interest  is  the  area  where  the  rapid  increase  in  end 
shortening  occurs.  We  identify  this  as  the  region  of 
instability. 

Since  the  magnitude  of  the  response  is  dependent 
on  the  mode  (m,n) , we  must  formulate  a criterion  which 
will  allow  us  to  choose  specific  m and  n for  a fixed- 
set  of  parameters;  ( ^ ) , ( , V'  , q.  We  propose  the 

following; 

The  correct  mode  is  represented  by  the  response 
curve  which  has  the  lowest  values  of  the  load  parameter 
in  the  region  of  instability.  This  criterion  is  com- 
patible with  a criterion  used  to  select  modes  in  the 
static  problem,  i.e.  minimizing  P^^^with  respect  to  the 
mode.  We  also  remark  that  this  mode  corresponds  to  the 
mode  for  which  the  amount  of  work  required  to  produce 
instability  is  a minimum. 

Thus,  selection  of  the  correct  mode  is  achieved  by 
comparing  the  response  curves  in  the  region  of  insta- 
bility and  choosing  that  curve  which  represents  a lower 
bound  on  the  load  parameter  in  the  region  of  instability 
(see  Figure  5) . This  may  appear  to  be  quite  a formidable 
task,  considering  the  number  of  reasonable  combinations 
of  m and  n available.  However,  one  can  reduce  this  number 


SOD  -XtajliS 
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e - dimensionless  end  shortening 


Figure  5.  Selection  of  Mode 
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considerably  by  the  following  procedure. 

Using  as  a starting  point  a number  of  circumfer- 
ential waves  slightly  higher  than  predicted  by  classical 
theory  for  an  equivalent  cylinder,  we  calculate  the  re- 
sponse for  some  range  of  m ( Lg  near  Ljp)  until  a lowest 
response  curve  is  found.  We  then  select  an  adjacent 
value  of  n,  and  repeat  the  process  for  m in  the  neigh- 
borhood of  the  previous  value  which  yielded  a minimum. 

We  continue  in  this  manner  until  the  lowest  curve  in  the 
region  of  instability  is  found.  It  should  be  noted  that 
the  load  in  the  region  of  instability  may  not  increase 
monotonically  for  all  values  of  m on  either  side  of  the 
lowest  curve;  for  in  a case  considered  herein,  another 
relative  low  curve  was  found  for  high  values  of  m.  Pur- 
suance of  this  relative  minimum  led  to  a ring-buckling 
type  failure,  with  the  number  of  longitudinal  waves  pre- 
dicted by  classical  theory. 

4 . 2 Effect  of  Loading  Rate 

We  shall  first  investigate  the  effect  produced  by 
various  loading  rates  on  a truncated  cone  having  the 
following  geometric  parameters 
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— = 1.171 

S| 


R9PI 

h 


1.75  X 10^ 


= .05  radians 


Since  we  have  neglected  the  kinetic  energy  of  the 
longitudinal  motion  in  our  analytical  formulation,  the 
loading  rate  is  restricted  by  the  following  consideration. 
The  time  to  reach  the  classical  buckling  load  should  be 
much  larger  than  the  time  required  for  a stress  wave  to 
travel  the  length  of  the  shell.  Symbolically, 


t 

CL. 


Sj-S, 


c 


where  c is  the  so-called  bar  velocity.  This  inequality 
leads  to  the  following  restriction  on  the  load  rate 

d^;  ( ^ _ 1 ) ctny 

In  our  particular  case. 


dt 


^ ^ .292 
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Thus,  computations  were  carried  out  for  the  following 
rates  of  loading 


dTf: 


.0050 


.0075  , .0100  , 


.0250  , .0500 


In  the  computations,  we  have  taken  U = .3  and  Uj,  = 1 x 10~^; 
the  latter  being  in-between  the  value  suggested  by  Loo  (19) 
and  Nash  (20) . The  results  are  presented  in  Figure  6, 
where  the  dimensionless  parameter  6 is  defined 


£ 


^ / ( 2 


ctn^^i* 

K 


) 


This  normalization  was  chosen  such  that  for  the  limiting 
case  of  a cylindrical  shell,  classical  buckling  occurs 
when  we  have  q = ?.  = 1.0. 

From  Figure  6 and  Table  1,  the  following  conclusions 
can  be  drawn.  Increasing  the  rate  of  loading  significantly 
increases  the  value  of  the  critical  load  in  the  region  of 
instability.  Also,  the  number  of  longitudinal  half-waves 
increases  with  increasing  rates  of  loading,  an  effect 
which  has  been  observed  in  experiments  on  cylindrical  shells. 
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It  is  also  apparent  that  one  cannot  assume  a fixed  ratio 
of  half-wave  lengths  for  the  solution  of  the  dynamic  prob- 
lem. The  tendency  toward  predominance  of  the  longitudinal 
waves  for  high  rates  of  compression  suggests  the  possi- 
bility that  for  rates  of  loading  in  which  the  longitudinal 
kinetic  energy  must  be  considered,  an  axisymmetric  formu- 
lation may  be  sufficient. 


TABLE  1 

RESULTS  OF  VARYING  THE  LOADING  RATE 


dq/d7S 

q 

^increase 

.0050 

.83-. 86 

— 

.81 

.0075 

.87-. 90 

5 

.90 

.0100 

.90-. 94 

8-10 

.99 

.0250 

1.02-1.10 

23-28 

1.23 

.0500 

1.15-1.28 

39-49 

1.90 
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4 . 3 Effect  of  Geometric  Parameters 

Although  of  interest,  the  case  of  holding  the  load 
rate  parameter  and  two  of  the  geometric  parameters  con- 
stant while  varying  the  third  does  not  represent  an  engi- 
neering type  situation.  Therefore,  this  procedure  has 
been  abandoned  in  favor  of  varying  the  geometry  in  a more 
physically  realizable  manner. 

Case  I .-  In  this  study  we  hold  the  total  height  (H) 

of  the  cone,  the  radius  of  the  cone  at  the  small  end,  and 
dp 

— constant,  while  increasing  the  conicity.  We  select  a 
cone  with  the  following  dimensions 

H = 6 inches 

Rj  = 1.75  inches 

h = .01  inches 

For  the  range  of  jT  considered,  the  values  of  the  parameters 
become 


sz/s, 

R^l  A 

d^dTt 

.050 

1.17 

175 

.0100 

.100 

1.35 

176 

.0101 

.200 

1.70 

179 

.0104 

.300 

2.06 

183 

.0110 

.400 

2.45 

190 

.0118 
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The  resulting  response  curves  are  shown  in  Figure  7, 

and  the  significant  results  tabulated  in  Table  2.  This 

study  indicates  a definite  increase  of  with  increasing 

t,  a fact  not  brought  out  in  a linear  analysis.  This  tend- 

Rd)i 

ency  (i.e.  for  constant  , the  load  parameter  increases 

h 

with  increasing  conicity)  has  support  in  the  experimental 
results  given  by  Seide  (2) . 

TABLE  2 

RESULTS  OF  VARYING  THE  GEOMETRY 
CASE  I 


“^CRlT. 

%increase 

4.r.  (1^^) 

in 

o 

• 

.90-. 94 

— 

205-214 

.99 

.10 

.92-. 96 

2 

208-217 

.94 

.20 

.97-1.02 

8-9 

212-223 

.97 

.30 

1.04-1.09 

15 

216-227 

1.25 

.40 

1.11-1.17 

23-25 

215-229 

1.22 

* For  a cone  fabricated  of  Mylar 


■^3  (1  ^ ) 2lTEh^  oos^t 
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1.4 


S - dimensionless  end  shortening 


Figure  7.  Effect  of  Varying  the  Geometry 

Case  I 
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Case  II  In  this  study  we  hold  the  total  height  (H) 

dp 

of  the  cone,  radius  of  the  cone  at  the  base,  and  -r-  constant, 

&& 

while  decreasing  the  conicity.  For  a cone  of  the  follow- 
ing dimensions 


H 

6.00  inches 

4.29  inches 

h 

0.01  inches 

and 

for  the  same 

range  of 

T as  used  in 

the  previous  study 

the 

values  of  the 

pertinent  parameters 

are 

t 

Sj/S, 

R<j9l  /h 

dq/dTS 

.40 

2.45 

190 

.0118 

.30 

1.76 

255 

.0110 

.20 

1.40 

314 

.0104 

.10 

1.16 

371 

.0101 

.05 

1.08 

400 

.0100 

The 

response  curves  are  shown  in  Figure 

8,  and  signifi- 

cant 

. results  in  Table  3. 

They  indicate 

a decided  de- 

0^  R^j. 

crease  in  q^^^^^with  increasing  ( — ) and  decreasing  con- 
icity.  Decreasing  q^^^_^with  increasing  ( — ) is  not 
surprising,  since  this  is  also  predicted  in  the  nonlinear 
analysis  of  cylindrical  shells. 
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n) 

(0 

O 

U 


w 

>= 

CN 


II 


1.4 


€.  - dimensionless  end  shortening 


Figure  8.  Effect  of  Varying  the  Geometry 

Case  II 
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TABLE  3 

RESULTS  OF  VARYING  THE  GEOMETRY 
CASE  II 


t 

R<pi/h 

rJ 

^CRlT. 

^decrease 

L^/Ls 

.40 

190 

1.11-1.17 

— 

1.22 

.30 

255 

.91-. 98 

18 

.91 

.20 

314 

.81-. 86 

27 

.91 

.10 

371 

.74-. 80 

32-33 

.88 

.05 

400 

.72-. 77 

34-35 

.73 

CHAPTER  V 


SUMMARY  AND  FURTHER  REMARKS 

The  equations  governing  the  nonlinear  response  of 
a truncated  conical  shell  to  time-dependent  axial  forces 
have  been  derived  by  employing  the  approximations  of  the 
nonlinear  theory  of  shallow  shells.  The  concept  of  in- 
itial imperfections  was  included  in  the  formulation. 

These  equations  can  be  easily  reduced  to  those  governing 
the  response  of  the  limiting  cases,  i.e.  circular  plate 
and  cylindrical  shell. 

A solution  form  for  conical  shells  of  sufficient 
length,  similar  to  that  used  in  previous  works  for  cylin- 
drical shells  and  indicated  in  experimental  results,  was 
assumed  and  the  Galerkin  method  used  to  obtain  the  differ- 
ential equations  which  govern  the  time-dependent  free 
parameters.  In  order  to  determine  the  amount  of  initial 
imperfection,  the  concept  of  an  "unevenness"  factor,  first 
used  by  Donnell  and  Wan,  was  introduced.  The  value  of 
this  factor  was  taken  near  those  used  successfully  in 
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previous  static  analyses  of  cylindrical  shells.  There- 
fore, to  obtain  reliable  quantitative  results,  experi- 
mental data  must  be  obtained  to  better  fix  this  quantity, 
and  to  determine  over  what  range  of  the  parameters  it 
can  be  considered  valid.  It  was  found  that  the  solution 
to  the  set  of  approximate  differential  equations  was  de- 
pendent on  three  geometric  and  one  physical  parameter. 
Because  the  assumed  solution  form  was  based  on  observa- 
tions of  cylindrical  shells  and  of  cones  of  conicity  less 
than  30°,  and  since  the  form  is  singular  at  ~ ~2’ 
approximate  solution  is  restricted  to  cones  of  moderate 
conicity. 

A criterion  was  proposed  to  determine  the  mode 
shape  of  the  buckled  cone.  It  was  compatible  with  that 
used  in  the  static  analysis,  and  yielded  a unique  mode 
form  in  all  cases  considered. 

In  the  investigation  of  the  effect  of  loading  rate, 
it  was  found  that  the  load  in  the  region  of  instability 
increased  significantly  with  increasing  load  rates,  as 
did  the  number  of  longitudinal  waves.  This  predominance 
of  longitudinal  waves  for  high-rates  of  loading  suggests 
the  possibility  that  for  the  case  in  which  the  longitu- 
dinal kinetic  energy  must  be  considered,  a nonlinear 
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axisymmetric  formulation  may  be  sufficient. 

From  the  studies  of  the  geometric  parameters,  the 
following  conclusions  can  be  drawn.  Increasing  the  con— 
icity  increases  the  load  parameter  in  the  region  of  in- 
stability, and,  thus,  the  case  of  a cylindrical  shell 
represents  the  lower  bound  for  thi^  critical  parameter. 
The  load  parameter  decreases  with  increasing  . The 
above  results  are  indicated  in  those  trends  obtained 
experimentally  by  Seide  (2).  However,  it  is  evident  that 
further  experimental  data  (static  as  well  as  dynamic) 
should  be  obtained,  guided  by  the  significant  parameters 
indicated  in  this  work. 

Although  the  numerical  studies  by  no  means  ex- 
hausted the  range  or  methods  of  varying  the  individual 
parameters,  it  is  felt  that  they  illustrated  the  impor- 
tant characteristics  of  the  response  and  the  trends  which 
are  evident  in  experimental  data,  but  which  are  not  in- 
dicated in  previous  theoretical  analyses. 


APPENDIX 


The  complete  Fortran  program  used  for  the  numerical 
solution  of  the  equations  which  govern  the  response  of  a 
truncated  conical  shell  to  a dynamic  axial  force  is  pre- 
sented on  the  following  pages, 

A liberal  number  of  comment  cards  have  been  intro- 
duced in  order  to  facilitate  correlation  with  the  main 
text.  It  is  also  noted  that  the  statements  have  been 
truncated  at  column  60  and  the  remaining  portion  listed 
on  the  following  line.  This  was  done  to  satisfy  margin 
requirements . 
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C FORTRAN  LANGUAGE  PROGRAM  FOR  SOLUTION  OF  THE 

C APPROXIMATE  GOVERNING  EQUATIONS  FOR  THE  NONLINEAR 

C RESPONSE  OF  A THIN  CONICAL  SHELL  TO  A DYNAMICALLY 

C APPLIED  AXIAL  FORCE 


CCTROL 

C CONICAL  SHELL  BUCKLING  CONTROL  PROGRAM 
DIMENSION  C(9),D(11),ZI5) 

COMMON  CtOf Z»T,SOS,RBYH,UOtAN,EM,EN,V,DTAUf F01.F02.CAY 

,El,E2,E3,ZO 

l,CTN 

C T=GAMMA,S0S=S2/Sl,RaYH=*RPHI/H,U0=IMPERFECTIUN  FACTOR 

C AN=POISSONS  RATIOt  EM=NO.  OF  HALF  WAVES  ALONG  THE 

C LENGTH,  EN=NO.  OF  WAVES  AROUND  THE  CIRCUM.,  A=ALPHA, 

C B=BETA,  CAY=K,  V=DQ/DTAU  WHERE  Q=Q-TILDA, 

C OTAU=DELTA  TAU 

18  READ  INPUT  TAPE  5 , 12 , T , SOS,RB YH , UO , AN , EM,EN, V, OTAU 
12  FORMAT  (7E10.3) 

CALL  GEODEV 
CONTINUE 
CALL  COEFT 
CONTINUE 
CALL  SOLVE 
CONTINUE 
GO  TO  18 
END 


CGEODV 

SUBROUTINE  GEODEV 
DIMENSION  C(9),D(11),Z(5) 

COMMON  C,D,Z,T,SOS,RBYH,UO,AN,EM,EN,V,DTAU,FOl,F02,CAY 

,El,E2,E3,ZO 

i,CTN 

CTN=COSF(T)/SINF(T) 

ZO=LOGF(SOS) 

CAY=SgRTF( 12.*{ l.-AN»*2) ) *RBYH* (C TN*»2 ) 

F01=UO*R8YH*( ( (SOS)**( l./£M)-l. )/( l.+( ( SOS ) ** ( 1 . /EM ) -1 

. )*EN/(3. 14* 

ISINFi T) ) ) )**2 
F02=.01*F01 
RETURN 
END 


CCOEFT 

SUBROUTINE  COEFT 

DIMENSION  C(9),D(11),Z{5),R(74),G(22) 

COMMON  C,D,Z,T,SOS,RBYH,UO,AN,EM,EN, V,DTAU,F01 ,F02,CAY 

,E1,E2,£3,Z0 

i,CTN 

C COMPUTATION  OF  REQUIRED  COEFFICIENTS 
A=3. 14159*EM/Z0 
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B=EN/SINF( T) 

AL=A»*2 

ALP=AL**2 

BE=B**2 

BET=BE**2 

Ull=(AL-2.+3.*(3.*AL+BE-4. > / ( AL+BE ) ) / ( ( AL+BE-4. } **2+16 

.*AL)*(-AU 

U12=( { AL-2. )*( l6.*BE-33. )-3.*{ ( 16 . *BE-9 . ) * ( BE-3 . *AL-4. 

)/{AL+BE)-6. 

i*(AL+BE-4.)))/(16.»( (AL+BE-4. )**2+16.*AL))*AL 
U21=(3.*AL-( AL-2. )*( 3.*AL+8E-4. ) / ( AL+BE ) ) / ( ( AL+8E-4 . ) * 

*2+16.*AL) *A 

U22=(-3.*AL*(16.*BE-33. )-{AL-2. )*{ ( 16. *BE-9 . ) * t BE-3. *A 

L-4. )/(AL+BE 

l)-6.*( AL+BE-4. } n/(16.*{ ( AL+BE-4. ) **2+16. »AL) ) *A 
U31=(6.*ALP-(2.*AL-1 . )*(3.*AL+1. ) )/( 16.*A*{ AL+1. )**2) 
U32=(-6.*A*ALP*BE-(2.*AL-1. )*(AL*(AL-1. )-(BE-l.)*(4.*A 

*AL-AL+1. ) ) ) 

i/( 32.*AL*{ AL+1. )**2) 

U33  = { (2.*AL-1. )*(3.*AL  + 1. )-6.*ALP) /( 128. *A* ( AL  + 1 . ) **2) 

*(4.*8E-3. ) 

U41=(3.*(3.*AL+1. )+2.*AL*(2.*AL-l. ) ) / ( 16 . * { AL+ 1 . )**2) 
U42=( 1.5*A*( AL*(AL-1. )-{BE-l. )*(4.*A*AL-AL+1. ) )-(2.*AL 

-1. )*ALP*BE) 

1/(16.*AL*( AL+1. )**2) 

U43=(-3.*{ 3.*AL+1. )-2.*AL*(2.*AL-l. ) ) / ( 128. * ( AL+ 1. ) **2 

)*(4.*B£-3. ) 

U5=AL/( 16.*(BE-l. ) ) 

U61=( (2.*AL-1. )*(3.*AL-BE+1. )/( AL+BE )-6 . *AL )/( 16 .*(( AL 

+BE-1. )**2+4 

1 .*AL ) ) *A 

U62=(2.*AL-1. )/( 32.* (AL+BE) ) *A 

U63=(6.*AL-(2.*AL-1. )*(3.*AL-BE+1. )/(AL+BE) )/(32.*( ( AL 

+BE-1. )**2+4 

l.*AL) )*A*(2.*BE-1. ) 

U71=(-3.*(3.*AL-8E+1. )/{AL+BE)-2.*(2.*AL-l. ) )/(16.*( (A 

L+BE-1. ) **2+ 

14.*AL) )*AL 

U72=-1.5*AL/( 16.* (AL+BE) ) 

U73=(1.5*(3.*AL-8E+1. )/(AL+BE)+2.*AL-l. )/( 16.*( (AL+BE- 

1. )**2+4.*AL 

1 )  )*AL*(2.*BE-1. ) 

U81=(-3.*(9.*AL-2. )*(8.*BE-1. )-9 . * ( 12 . *AL* ( 24. *BE-27. ) 

-(9.*AL+B£-4 

1. )*( 10.*8E+18.*AL-9. ) )/(9.* AL+BE) )/( 16 . * ( ( 9 . *AL+BE-4 . ) 

**2+144. *AL) 

2) *AL 

U91=(27.*AL*(8.*BE-l . )-(9.*AL~2. ) * ( 12. *AL* ( 24. *BE-2 7 . ) 

-(9.*AL+BE-4 

1. )*( 10.*BE+18.*AL-9. ))/(9.*AL+8E) )/( 16.*( ( 9.*AL+BE-4. ) 

**2+144. *AL) 

2)*A 

U101=(-( AL-2. )*(24.*8E*( AL+9.*BE-4. )-( 3.*8E+11 .*AL) )-3 

.*( 12.*AL*(8 
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)-(AL+9.*8E-4. )»(6.*8E-2.«AL-3. ) ) )/ ( l6.*(AL+9.* 

BE)*( (AL+9.* 

2BE-4. )**2+16.*AL) ) *AL 

Uill=l 3.*AL*(24.*BE»(AL+9.*BE-4. ) - ( 3 .*BE+l 1 . *AU )-( AL- 

2.  )*( 12.*AL» 

i(8.*BE-l. )-(AL+9.*8E-4. )*(6.*BE-2.*AL-3. ) ) ) / ( 16. *( AL+9 

.*BE)*( (AL+9 

2.»BE-4. )**2+16.*AL ) )*A 

U121={-(9.*AL-2. )-3.»(2.*(9.*AL+9.*BE-4. )+ { 27. *AL-9 . »B 

E-4. }/(AL+BE 

1)  ) )/( 16.*( (9.*AL  + 9.»BE-4. ) **2  + 144 . *AL ) )»AL 
U13i=(9.*AL-l./3.*(9.»AL-2. J * (2. * (9. *AL+9. *8E-4. )+(27. 

*AL~9 . *BE+4. 

1 )/ ( AL+BE) ) ) /( 16.*( (9.*AL+9.*BE-4. )**2+144.*AL) J *A 
U141=(-1./16.*(8.*AL-1. ) *(48.*(8E-AL )+(4.*BE-24.*AL-3. 

)/AL)+24.*AL 

i*BE)/(64.*( (4.»AL-1. )**2+l6.*AL) )*A 
U151=( 3./8.*{48.*(BE-AU+(4.*8E-24.*AL-3. )/AL) +4.*(8.* 

AL-1.  )*BE)/( 

164. *( (4.»AL-1. )**2+16.*AL) )*AL 
U161=-AL/( 16.*{ 16.*8E-l. ) ) 

U171=(-.25*(8.*AL-1. )*(2.*(8.*ALP-2.*AL*BE+BET}+(8.*AL 

-3.*BE+1.) )/ 

l(4.»AL+BE)-6.*AL»BE)/{ 16.* ( (4.*AL+BE-1. )**2+4.*AL) ) *A 
U181=( 1.5*(2.*I8.*ALP-2.*AL*BE+BET)+(8.*AL-3.*BE+1. ) )/ 

(4.*AL+BE)-{ 

18.*AL-1. )*BE)/{ 16. *( (4.*AL+8E-1. )**2+4.*AL) ) *AL 
U 191= (-(2.* AL-1. )* (4.*BE-3.*AL-1 . ) / ( AL+4 . *BE ) -6. *AL ) / ( 

8.*16.*( (AL+ 

14.*BE-1. )**2+4.*AL) ) *A* ( 4. *BE-1 . ) 

U201=( 3.*(4.*BE-3.*AL-1. ) / ( AL+4 . *8E ) -2 . * ( 2 .*AL- I . ) )/(8 

.*16.*( (AL+4 

l.*BE-l. )**2+4.*AL) ) *AL*(4.*b£-l. ) 

U211=(8.*AL-1. )*A/(32.**2*(AL+BE) ) 
U221=-6.*AL/(32.**2*(AL+BE) ) 

(311  = ( ( 8E-2.  ) + (3.*AL+BE-4.  )/(  AL+BE)  ) / ( ( AL+BE-4.  )**2+16. 

*AL)*(-AL) 

Q12=( ( BE-2. )»( 16.*BE>33. )-( ( 16. *BE-9. ) * ( BE-3. * AL-4 . )/( 

AL+BE)-6.*(A 

lL+BE-4. ) ) )/( 16.*( (AL+BE-4. )»*2+16.*AL) ) »AL 
Q21=(AL-(BE-2. )*(3.* AL+BE-4. )/(AL+8E) )/( ( AL+BE-4. ) **2+ 

16.*AL)*A 

Q22=(AL*( 16.*BE-33. )+(BE-2. )*( ( 16. *BE-9. ) * ( BE-3. *AL-4. 

)/(AL+BE)-6. 

I* (AL+BE-4. ) ) )/( 16.*( (AL+BE-4. )**2+16.*AL) )*(-A) 
Q31=(2.*ALP+3.*AL+1. )/( 16.*A*(AL+1. )**2) 

Q32=( AL*(AL-1. )-(BE-l. ) * (4. *A*AL-AL+ 1. ) -2 . *A*ALP*BE ) / ( 

32.*AL*( AL+1 

1. )**2) 

Q33=(2.*ALP+3.*AL+1. )/( 128.«A*(AL+1. ) **2 ) * ( 3 .-4 . *BE ) 
Q41=l./( 16.*( AL+1. ) ) 

Q42=( (AL*(AL-1. )-(BE-l. ) * ( 4 . * A*AL-AL+ 1 . ) ) +2 . *A*AL*BE ) / 

(32.*A*( AL+1 


1. ) **2) 
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Q43=  (3.-4.*BE)/(128.»(AL+l.n 

Q5=( l.-2.*BE)*AL/( 16.»(BE~1.  ) ) 

Q61  = { t 2.*BE-l. )*(3.*AL-BE-l. j/(AL+BE)-2.*AL)/( 16.*C (AL 

+8E-1. )**2+4 

1. *AL) )*A 

Q62=(2.*8E-l. )*A/(32.*(AL+BE) ) 

Q63=(-( 2.*B£-1. ) •{ 3.*AL-BE+1. ) / ( AL+BE ) +2. *AL ) / ( 32. » ( (A 

L+BE-1.)**2+ 

14.*AL) )»A»{2.*BE-1.) 

Q71=( (3.»AL-BE+1. )/{AL+B£)+2.*(2.*BE-l. ) )/ ( i6.»( (AL+BE 

-1.  )*»2  + 4.»A 

ID  )*(-AD 

Q72=-A/(32.»(AL+BE) ) 

Q73=( (3.*AL-8£+l. )/( AL+BE)+2.*(2.*BE-1. ))/ (32. *( (AL+BE 

-I . ) »*2+4. »A 

ID )*AL*(2.*BE-1. I 

Q81=(3.*(BE-2. )*(8.*8E-1. )+3.»( l2.*AL»(24.*BE-27. )-(9. 

•AL+BE-4. )*( 

ilO.»BE+l8.*AL-9.) )/(9.»AL+BE) ) / ( 16 . * ( ( 9 . *AL+BE-4 . ) *»2+ 

144. *AD  )*(- 

2AL) 

a91=(9.*AL*(8.»BE-l. )-(BE-2. )*(12.*AL*(24.*8E-27. )-(9. 

•AL+BE-4. ) »( 

110.*BE+18.*AL-9. ) )/(9.*AL+BE) )/( 16.*( (9.»AL+BE-4.)**2+ 

144. *AD  )*A 

Q101=( (9.*BE-2. ) •(24.*8E*( AL+9.*BE-4. )-(3.*BE+ll.*AL ) ) 

+( 12.*AL*(8. 

l*8E-l. )-(AL+9.*BE-4. ) • ( 6 . »B£-2 . * AL-3 . ) ) ) / ( 16 . * ( AL+9 . *B 

E)*( (AL+9.*B 

2E-4.  )*»2  + 16.*AD  )*(-AD 

Qll 1=( AL»(24.*BE*( AL+9.*BE-4. )- ( 3 . *BE+ 1 1 . *AD )-(9.*BE- 

2.  )•( 12.*AL* 

L(8.»BE-i. )-( AL+9.*BE-4. ) * ( 6 . *BE-2. »AL-3 . ) J )/ ( 16.* (AL+9 

.*BE)*( (AL+9 

2. *BE-4. )**2+16.*AD )*A 

Q121=( (9.*BE-2. )+(2.*(9.*AL+9.*BE-4. )+ ( 27. *AL-9. *BE+4. 

)/(AL  + BD)  )/ 

1(16.*( (9.*AL+9.*BE-4. )**2+144.*AD )»(-AL) 

Q131=( 3.*AL-l./3.*(9.*BE-2. ) *(2.*(9.*AL+9.*BE-4. )+(27. 

*AL-9.*BE+4. 

1 ) / ( AL+BE ) ) )/ ( 16.*( (9.*AL+9.»BE-4. ) **2+144.* AD ) *A 
Q141=( l./16.*(4a.* (BE-AD+(4.*BE-24.*AL-3. )/AL)+8,»AL* 

BE)/(64.*( (4 

i.*AL-l. )**2+16.*AL ) )*A 

Q151=(  l./8.*(48.*(BE-AD  + (4.*BE-24.*AL-3.  )/AD-4.*BE)/ 

(64. *( (4.*AL 

1-1  . )**2+16.*AL ) ( *AL 
Q161=(8.*BE-l. )*AL/( 16.* ( 16.*8E-1. ) ) 

U17l=(-.25*(2.*BE-1. )*(2.*( 8.*ALP-2.*AL•BE+BET)+8.*AL- 

3. *BE+l.  )/(4 

l.*AL+8E)-2.*AL*8E)/( 16.* ( ( 4. *AL+BE-1 . ) **2+1 6. *AD ) *A 
Q181=( .5*(2.*( 8.*ALP-2.*AL*BE+BET)+8.*AL-3.*BE+1. )/(4. 

*AL+BD-BE*( 

12.*BE-1. ) J/( 16.*( (4.*AL+BE-1. )**2+16.*AD )*AL 
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Qi91=(-{8.*BE-1. )*(4.»BE-3.*AL-1. ) / ( AL+4.*BE )-2.*AU / ( 

128.« ( (AL+4. 

l*BE-l. )**2+4.*AL ) ) »A*(4.*BE-1. ) 

Q201=( (4.*BE-3.*AL-i. )/{AL+4.*BE)-2.*(8.*BE-l. ) )/( 128. 

*( (AL+4.*8E- 

11. )**2+4.*AL) )*AL*(4.*BE-1.) 

Q211=(8.*BE-1. )*A/(32.**2*(AL+BE)) 

Q221=-2.*AL/(32.**2*( AL+BE) J 

Vll=( ( 3.*AL+BE-4. )/(AL+BE)-l. )/( (AL+BE-4. ) **2+16.*AU* 

AL*B 

V12=( 16.*BE-33.+{ (16. »BE-9. )*(BE-3.*AL-4. )/{ AL+BE )~6.» 

(AL+BE-4. ) ) ) 

l/i (AL+BE-4. )**2+16.*AL)*AL*B/16. 

V21=(AL+(3.*AL+BE-4. )/(AL+BE))/( (AL+BE-4. ) **2+16.* AL )* 

(-A)*b 

V22=( AL*( 16.*BE-33. )-( 16 . *BE-9. )*( BE-3. *AL-4. )/( AL+BE) 

+6 . * ( AL+BE-4 

!.))/( (AL+BE-4. )**2+i6.*AL)»A*B/I6. 

V5=AL*B/(16.*(BE-1.) ) 

V61=(4.*AL+( 3.*AL-BE+1. )/(AL+BE) )/( 16.* ( ( AL+BE-1 . ) **2+ 

4.*AL) )*(-A) 

i*B 

V62=-A*B/( 32.* (AL+BE ) ) 

V63=(4.*AL+(3.*AL-8E+1. )/( AL+BE) )/(32.*( ( AL+BE- 1 .) **2+ 

4.*AL  ) )*A*B* 

1(2.*8E-1. ) 

V71=(1.-(3.*AL-BE+1. )/(AL+BE) )/(8.*( ( AL+BE-i . ) **2+4. *A 

L) )*AL*8 

V72=-AL*B/( 16.* (AL+BE) ) 

V73=( ( 3.*AL-BE-1. )/( AL+BE )-l. )/( 16.* ( (AL+BE-1. ) **2+4.* 

AL) )*AL*B*(2 

1. *BE-l. ) 

V81=(-3.*(8.*BE-1. )+3.*( 12.*AL*(24.*BE-27. )-(9.*AL+BE- 

4. )*( 10.*BE+ 

118. *AL-9. ) )/ (9.* AL+BE) )/ ( 16.* ( (9.* AL+BE-4. ) **2+144. »AL 

) ) *AL*B 

V91=(9.*AL*(8.*BE-1. )+ ( 12. *AL* ( 24.*BE-27. )- ( 9. *AL+BE-4 

. )*( 10.*BE+1 

18.*AL-9. ) )/(9.*AL+B£) )/( 16.*( (9. »AL+BE-4. ) **2+144. *AL) 

)*(-A)*B 

V101=( ( 12.*AL*(8.*BE-1. )- ( AL  + 9.*BE-4. ) * ( 6. *8E-2. *AL-3. 

) )-(24.*BE*( 

lAL+9.*BE-4. )-(3.*BE+ll.*AL) ) ) / ( 16. • ( AL+9.*BE ) * ( ( AL+9.* 

BE-4. ) **2+16 

2. *AL) )*3.*AL*B 

Vlll=( AL*(24.*BE*( AL+9.*BE-4. )- ( 3. *BE+ 1 1 .*AL ) ) + ( 12  - *AL 

*(8.*BE-1.)- 

1 (AL+9.* BE-4. )*(6.*BE-2.*AL-3. ) ) ) / ( 16 . * ( AL+9. *BE ) * ( (AL+ 

9.*BE-4. )**2 

2+l6.*AL) )*3.*{-A)*B 

V121=(2.*(9.*AL+9.*BE-4. )+( 27 . *AL-9. *B£+4. )/( AL+BE ) -1 . 

)/(16.*( (9.* 

lAL+9.*B£-4. ) **2+144. *AL) )*3.*AL*B 


61 


Vl3i=(2.*(9.*AL+9. *BE-4. ) + ( 27. »AL-9 . *BE+4. )/ ( AL+BE) +9. 

•AL)/(16.*( ( 

19.*AL+9.*BE-4. ) **2+144. *AL) )*(-A)*B 
V161=-AL*B/(8.*(16.*8E-1.)  ) 

V171=( .25*{2.*l8.*ALP-2.*AL*BE+8ET)+(8.*AL-3.*BE+l. ) )/ 

(4.*AL+BE)-4 

i.*AL*8£)/( 16. *( t4.*AL+BE-l. )**2+16.*AL) )*A*8 
V181=(BE+(2.*(8.*ALP-2.*AL*BE+8ET)+(8.*AL-3.*BE+1. ) )/( 

4.*AL+BE) )/( 

116.  *{ (4.*AL  + 8E-1. )»*2+16.*AL) )*AL*B 
V191=l .25*(4.*BE-3.*AL-1. )/(AL+4.*BE)-AL)/( 16.*{ (AL+4. 

*BE-1 . ) **2+4 

l.*AL) )*A*B»(4.*BE-1. ) 

V201=( l.+(4.*BE-3.*AL-l. )/(AL+4.*BE) )/(32.*( (AL+4.*BE- 

1.  )**2+4.*AL 

1)  )*AL*B*(4.*BE-1. ) 

V211=-A*B/(512.*(AL+B£) ) 

V221=-AL*8/{ 128.*( AL+BE) ) 

R( 1 ) =.25*( (BE-1 . )»(2.*U41-U71)+AL*(2.*G41-Q71 )+A*(2.* 

(U31+Q31)-(U 

161+Q61) )+2.*A*B*V61) 

R(2  ) =.25* ( (BE-l. )*(2.*U42+2. *U5-U72 ) +AL* ( 2 . *Q42+2. *G5 

-Q72)+A*(2.* 

1(U32+Q32)-(U62+Q62) ) +2 . * A*B* V62 ) 

R(3  ) =.25* ( ( BE-l. )*(2.*U43+2.*U5-U73  J + AL*(2.*Q43+2.*Q5 

-Q73)+A*(2.* 

1(U33+Q33)-(U63+Q63) )+2.*A*B*V63) 

R(4)  =3./16.*((4.*BE-l. ) *U1 1+4. *AL*Q1 1-2 .* A* ( U2 1+Q2 1 )- 

2. *Ull-16./6 

l.*A*8*V21) 

R{ 5 ) = .0625*(  ( 4.*BE-1. ) * ( 3 .*U101-3.*U12  + U81-U12U + 4. *A 

L*{3.»Q81-3. 

1*Q12+Q101-Q12U+2.*A*(3.*(U22+Q22)+(U131+Q131)-3.*(U91 

+Q91)-(U111+ 

2Q111) )+2.*( 3.*U12-U81)+8.*A*B*( V22-V91-V111+V131) ) 

R( 6 ) =.2  5*( ( BE-i. )* {U61-2.*U31)  + AL*(Q61-2.*Q31 )+A* (2.* 

(U41+Q41)-(U 

171+Q71 ) )+2.*A*B*V71) 

R( 7)  =.25*< { BE-l. ) *( U62-2.*U32)+AL*{ Q62-2.*Q32 )+A* (2.* 

(U42+Q42)-2. 

1*(U5+Q5)-(U72+Q72) ) +2 . * A*B* ( 2. *V5+V72 ) ) 

R{ 8)  =.25* { (BE-l. )*(U63-2.*U33)+AL*( Q63-2.*Q33  )+A*( 2.* 

(U43+Q43)-2. 

1*(U5+Q5)-(U72+Q72) ) +2 . *A*B* ( 2 . * V5+ V73 ) ) 

R(9)  =1./I6.*l (4.*BE-l. )*U2i-l2.*AL*Q21-6.*A*(Ull+QllJ 

-8.*A*B*Vll- 

12.*U21) 

R(l0)=l./16.*{ (4.*BE-1. )*(U91-U22+Ulll-U131)+4.*AL*(3. 

*Q22+3.*Q91- 

IQll 1-Q121 )+2.*A*(3.*{ U12+Q12)+3. *(U81+Q81 )-( UlOl+QlOll 

+(U121+U121) 

2) +8.*A*B*( V12+V81-V101-V121 )+2.* (U22-U91) ) 

R( 11  ) = .25*( (BE-l. )*(U71-2.*U41)  + AL*(Q71-2.*Q41 )+A*(2.* 

(U31+G31 )-(U 


62 


161+G61) )+2.*A*B«V61) 

R(12)=.25*( ( BE-1. )*(U72-2.*U42)+AL*(Q72-2.*Q42)+A*(2.» 

(U32+Q32)-(U 

162+C62 ) ) +2.*A*B*V62I 

R( 13)=.25*( ( BE-1. )* (U73-2.*U43+2.*Ui51-U181 ) +AL» (Q73-2 

.*Q43+2.*Q15 

11-Q181 )+A*(2.*(U33+Q33)+2.*{U141+Q141)-(U63+Q63)-(U17i 

+Q171) )+2.*A 

2*8*(V63+V171) ) 

R(14)=l./16.*( (4.*BE-1. )*Ull+12.*AL*Qil+6.*(U21+Q21 )-2 

.•Ull+8.*A*B 

l*V21) 

R(  15)  = 1./16.*(  (4.*BE-1.  )*{Ui2-U10l-2.*U8H-Ul21)  + 4.*AL* 

(3.*Q12-Q101 

1) +2.*A*( 3.*{U22+Q22)-(Ulll+Qlll) ) +2 . * ( 2. *U81-U12 ) +8. »A 

*B*( V22+V111 

2)  ) 

R{ 16)=.25*( (BE-l. )*i2.*U31-U61)+AL*{2.*Q3i-Q6l )+A*(2.* 

(U41+Q41 )-(U 

171+071) )+2,*A*B»V71) 

R(  17)  = .25*( (BE-1. )*(2.*U32-U62)+AL*(2.*Q32-Q62)  + A*(2.* 

(U42+042)-(U 

172+Q72) )+2.»A*B*V72) 

R( 18)=.25*( ( BE-1. )*( 2.*U33-U63-2.*U141+U171 )+AL»(2.*C3 

3-063-2. *014 

ll+Ql71)+A*(2.*(U43+043)-{U73+Q73)+2.*(U151+Q151)-(U181 

+0181) )+2.*A 

2*B*(V73+V181) ) 

R(19)=1./16.*((4.*BE-1. )*U21+12.*AL*021-6.*A*IU11+011) 

-8.*A*B*V11- 

12.*U21) 

R(20)=l./16.*( (4.* BE-1. )*(U22-Ulll-2.*U91+2.»U131)+4.» 

AL*(3. *022-0 

1111)+2.*A*((U101+Q101)-3.*IU12+012))+8.*A*B*(V101-V12) 

+2.*(2.*U91- 

2U22  ) ) 

R( 21)=.25*( (BE-1. )*{U181-2.*U151)+AL*{Q181-2.*Q151 )+A* 

(2.*(U141+Q1 

141 )-(U17l+0171) )+2.*A*B*V17l) 

R(22)  = l./16.*(  (4.* BE-1.  ) * ( U81-II121 ) +4 . *AL*  ( 3 . *08 1-0121 

)+2.*A*( 3.*( 

1U9 1+091 )-( U1 31+0131 ) )+8.*A*B* (V91-V131 )-2.*U81 ) 
R(23)=.25*((BE-1.)*(2.*U141-U171 )+AL*{ 2. *0141-0171 )+A* 

(2.*(U15l+01 

151 )- (1)181+0181 ) )+2.*A*B*V18l ) 

R(24)=1./16.*( (4.* BE-1. )*(U91-U131)+4.*AL*(3.*09l-0131 

)+2.*A*( (U12 

11+Q121)-3.*(U81+Q81 ) )+8.*A*B*(V12l-V81)-2.*U91 ) 

R( 25)=.0625*( (4.*BE-1. ) *U71-4 .*AL* ( 2 . *04 1-071 ) -2 .* A* (2 

.*(U31  + 03D- 

1 (U6 1+061) )-8.*A*B*V61+2.*U41) 

R(26)=.0625*( (4.*BE-i. ) * ( 2 . *U5-U72 ) +4. *AL* ( 2 . *042-072 ) 

+2.*A»(2.*(U 

132+032 )-(U62+062) ) +8 . *A*B*V62-2. *U42 ) 
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R(  27)  = .0625*(  U.*BE-l.  ) * { 2. *U5-U73 ) +A. *AL*  ( 2 .*Q43-Q73 ) 

+2.*A*(2.*(U 

133+Q33)-(U63+Q63) )+8.*A*B*V63-2.*U43 J 
R(28)=.25*( {BE-1.)*U12+AL*Q12-A*(U22+Q22)-2.*A*B*V22) 
R(29)=.25*( (BE-1. )*Ull+AL*Qil-A*{U21+Q2l )-2.*A*B»V21) 
R(30)=.125»( {4.*BE-i. >*U61+2.*U31) 

R(3U=. 125*( (4.*BE«1. )*U62-2.»A* (U5+Q5 )+8.*A*B»V5+2.*U 

32) 

R(  32)  = .0625*( (4.»BE-1. )*(2.*U63-U171)+4.*AL*{2.*G141-Q 

171)+2.*A*( ( 

IU181+Q181)-2.*(U5+Q5+U151  + Q151) J +8 . *A*B* ( 2 . » V5-V181 ) +2 

.*(2.»U33-Ul 

241 ) ) 

R{ 33)=.25*( lBE-1. )*(U91-U22)+AL*{Q91-Q22)+A*(U12+Q12+U 

8l+Q8i)+2.*A 

1*B*{V12+V81) ) 

R(34)=.25*( (BE-1. )*U21+AL*Q21-A*(U11+QII)-2.»A*B»V1U 
R(35)=.125»< {4.»BE-1. )*U71+2.»U41) 

R(36)=.125*( (4.»BE'1. )*(U72-U5)-4.*AL*Q3+2.*U42) 

R{ 37) =.0625* ( (4.*BE-1. )* ( 2.*U73-2.*U5-U181 ) +4.*AL*( 2.* 

Q151-2.*05-Q 

1181)+2.»A*(2.*(U141+Q141 )-(U171+Q171 ) ) +8.*A*B*V171+2.* 

(2.*U43-U151 

2)  ) 

R{ 38)=.25*( ( BE-1. ) * ( U12-U81 ) +AL* ( Q12-Q8i ) +A* ( U22+Q22+U 

91+Q9i)+2.*A 

1*B*(V22+V91) ) 

R(39)=.25*( (BE-1. )*Ull+AL*Qil+A*(U21+Q21)+2.*A*B»V2l) 
R(40)=.25*( (4.*BE-1. )*U31+U61 ) 

R(41)=.25*( (4.*BE-1. )*U32+2.*A*(U5+Q5)+U62) 

R(42)  = .0625*(  (4.*BE-1.  )*(  4. *U33+2.  *1)191-2.  *U  14 1-U2  1 1 ) + 

4.*AL*(2.*Q1 

i71-2.*Q141-Q211)  *^2.*A*  (4.*(U5+Q5  )-2.*(U181+Q181)+2.»  (U 

151  + Q15D-2. 

2*(U161+Q161)+(U221+0221 ) ) +8. *A*B* ( 2 . *V 161-V22 1 ) +2. • ( 2. 

•U63-U171) ) 

R(43)=.25*( (BE-l. )*( U22-U91-U111+U131)+AL*(G22-Q91-Q11 

1+G131)+A*(U 

1101+Q101+U121+Q121-U12-Q12-U81-Q81)+2.*A*6*( V12+V81+VI 

01+V121) ) 

R(44)=.25*( (BE-1.)*U21+AL*G21-A*(U11+Q11)+2.*A*B*V11) 
R(45)=.25*( (4.*BE-1. )*U41+4.*AL*Q41+G71) 

R(46)=.25*( (4.»BE-1. ) *U42+4. *AL*Q42^ ( U72-U5) ) 

R( 47) =.062 5* ( (4.*BE-1. )*(4.»U43+2.*U201-2.*U151-2.*U16 

1-U221)+4.*A 

1L*( 4.*Q43+2.*Q181-2.*G151-2.*G161-Q221 ) +2.*A*(2.»(U171 

+Q171)-2.*(U 

2141 *0141 )-(U211+Q2 11 ) ) +8 . *A*B*V2 1 1+2 . * ( 2 . *U73-2 . *U5-U1 

81)  ) 

R(48)=.25*{ (BE-l. )*(U81-U12+U10i-U121)+A*( (U111+Q111)+ 

(U131+Q131 )- 

i(U22+Q22)-{U91+G91 ) ) +2 . * A*B* ( V22+V91+V1 1 1+Vl 31 ) +AL* ( Q8 

1-G12+Q101-Q 


2121)) 
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R(49)  = .25*(  (BE-l.  )*UiH-AL*Qli+A*(U2l  + Q2l)-2.*A*8*V2l) 
R(50)=.0625*( (4.*BE-1. )*U61-4.*AL*(2.*Q31-Q6l)-2.*A*(2 

.*  (U4H-Q41  )- 

1 (U71+Q71 ) )-8.»A*B*V71+2.»U3l) 

R{ 51 )=.0625* I ( A.*BE-1. ) *U62-4.*AL* { 2 .*Q32-Q62 ) -2 .*A* ( 2 

.*(U42+QA2)- 

UU72+Q72) )-8.»A*B*V72+2.»U32) 

R(52)=.0625*( (4.*BE-1. )»(2.»U171-U63)+4.*AL*(2.*Q33-Q6 

3)+2.*A*(2.* 

1(U43+Q43)-{U73+Q73) )+8.*A*B*V73+2.«(2.*U141-U33) ) 

R( 53)=.25*( (BE-1. )*U91+AL*Q91-A* (U81+Q81 )-2.*A*B*V81) 
R(54)=.0625»{ (4.*BE-1. ) *U71-4. »AL* ( 2 .»Q4l-Q7 I ) -2.* A* ( ( 

U61+Q6D-2.* 

UU31+Q31) )+8.*A*B*V6l+2.*U41) 

R(55)=.0625*( (4.*BE-1. )*U72-4.*AL» ( 2 .*Q42-Q72 ) -2 . • A» ( ( 

U62+Q62)-2.* 

1(U32+Q32))+8.*A*B*V62+2.*U42) 

R(56)=.0625»( (4.*BE-1. )*(2.*U181-U73)+4.*AL»(2.»Q43+Q7 

3)+2.*A*( (U6 

13+Q63)-2.*(U33+Q33) )-8.»A*B*V63+2.*(2.*U151-U43) ) 

R( 57)=.25*{ (BE-1. ) *U81+AL*Q8 i+A* ( U91 +Q91 ) +2. *A*B*V91 ) 
R(58»=.125*( (4.*BE-1. )*U31-4.*AL*(Q61-Q31)-2.*A*( (U71+ 

Q71)-(U41+Q4 

11)  )+U61) 

R(59)  = .125*( (4.*BE-1. )*U32-4.*AL*(Q62-Q32)-2.*A»{  (U72+ 

Q72)-(U42+Q4 

12)  )+U62) 

R(60)=.0625*( (4.*BE-1. )*(4.*U141+2.*U211-2.*U33-U191)+ 

4.*AL*(2.*Q6 

13-2.*Q33-Q191 )+2.*A*(2.*(U73+Q73)-2.*(U43+Q43)-(U201+Q 

201) )+8.*A*B 

2*V201+2.*(2.*U171-U63) ) 

R(61)=.25*( (BE-i. ) *(U91-U131)+AL*(Q91-Q131 )+A*( (U121+Q 

121)-(U81+Q8 

11)  )+2.*A»B*( V81+V121) ) 

R(62)=.125*{ (4.*BE-1. )*U41-4.*AL»(Q71-Q41)-2.*A»( (U31+ 

Q31)-(U61+Q6 

11)  )+U71) 

R( 63)=.125*( (4-*BE-l. )*(2.*U5-U42)+4.*AL»(Q72-Q42)+2.* 

A*( (U32+Q32) 

1-(U62+Q62) )-U72) 

R(64)=.0625*( (4.*BE-1. ) * ( 4. *U5+2 . *U221-2. *U43-U201 ) *4. 

*AL*(2.*Q73- 

12.*Q43-Q201)+2.*A*(2.*(U33+Q33)+(U191+Q191)-2.*(U63+Q6 

3) )-8.*A*8*V 

2191+2.*(2.*U181-U73) ) 

R(65)=.25*( (BE-1. )*(U121-U81)+AL*(Q121-Q81)+A*( (U131+Q 

131)-(U91+Q9 

11)  )+2.»A»B*(V91+V13i ) ) 

R(66)=.0625*( (4.*BE-1. )*U171-4.*AL*(2.*U141-U171)-2.*A 

*(2.»(U151+Q 

1151)-(U181-»-Q181)  )-8.*A*8*U181+2.*U141) 

R(67)  = .0625*(  (4.*BE-1.  ) »U18l-4.  *AL»  ( 2 .*(il51-Q181  )-2.»A 

*{ (U171+Q171 
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1 )-2 .*(U141+Q141) ) ♦8.*A*B*V171+2.»U151) 

R( 68) =.0625* ( (4.»BE-1. )* ( U2 i 1+2 . *U14 1 ) -4 . *AL*( 2 .*Q171- 

2. *0141-0211 

1 )-2.  »A* ( 2. * I U18 1 + 0181) -2.* (U 151+0151 )-( U22 1 + 0221 )) -8. * 

A*B*V221+2.* 

2U171) 

R(69)=.0625*( (4.*8E-1. )* (U221+2.*U151)-4.*AL*(2.*C181- 

2.*0151-0221 

I) -2.*A*{2.*(U141+0141)+(U211+0211)-2.*(U171+Q171) )+8.* 

A*B*V211+2.* 

2U181) 

R(70)=.125*{ (4.*BE-1. ) *U41-4 . *AL* ( 071-Q41 )-2.*A* (I U61+ 

061) -{U31+03 

II)  )+U71) 

R(  71)  = . 125* ( (4.*BE-1. )*U42-4.*AL»(072-Q42)-2.*A* ( (U62+ 

062) -(U32+03 

12) )-(2.*U5-U72) ) 

R{ 72)=. 0625* ( (4.*BE-1. ) * ( 2. *U161-2 .*U43-U201 ) +4. *AL* (2 

.*073-2. *043 

1-0201 )+2.*A* (2.* (U63+063 )-2 . * (U33+Q33 ) - ( U191+Q191 ) )+8. 

*A*B*V191+2. 

2*(2.*U5-U73) ) 

R{ 73)=.25*( ( BE-1. )*(U12-U101 )+AL*{012-Q101 )+A*( (Ulll+Q 

lll)-lU22+02 

12) ) +2.*A»B*( V22+V111 ) ) 

R( 74)=.25*( (BE-1. )*U11+AL*Q11-A*(U21+Q21 )+2.*A*B*V21) 
G( 1 )=AL/(AL  + 1. ) 

G(2  ) = A/{ AL+1. ) 

G(3)=3.*AL/(4.*ALP+5.*AL+1. ) 

G(4)=A*(2.*AL-1. )/(4.*ALP+5.*AL+l. ) 
G(5)=5.*AL/(36.*ALP+13.»AL+1. ) 

G(6  ) = A* (6.*AL-1. )/(36.*ALP+13.*AL  + l. ) 

G(7)=AL/(AL+1. ) 

G( 8)=3.*A*AL/(4.*ALP+5.*AL+i. ) 

G(9)=AL*(2.*AL-1. ) / ( 4. *ALP+5 .*AL+ 1 . ) 

G( 10)=6.*A*AL/( (4.*AL+1. )*(9.*AL+1. ) ) * ( AL- I . ) / ( AL+ 1 . ) 
G( 11)=AL*( ll.*AL-l. )/{ (AL+1. ) * ( 4. * AL +1 . ) * ( 9 . *AL+1 . ) ) 

G( 12)=3.*A*AL/( (4.*AL+1. )*(9.*AL+1. ) )*(8.*AL-3. )/( 16.* 

AL+1. ) 

G( 13)=AL*(26.*AL-1. ) /( { 4 . * AL+ 1 . ) * ( 9. *AL+ 1 . ) * ( 16. *AL+ 1 . 

) ) 

G( 14)=ALP+BET+2.*AL*BE+2.*AL-2.*BE+1. 

G( 15)=AL*(BE-AL-1. ) 

G( 16)=2.*AL*(AL+BE+1. )/(AL+l. ) 

G(17)=AL+BE+1. 

G( 18)=8E-AL-1. 

G( 19)=AL/(4.*AL+9. ) 

G( 20)=. 125*(48.*ALP+48.*BET+32.*AL*BE-24.*BE+40.*AL+9. 

) 

G(21)=. 125*AL*( 12.*BE-12.*AL-9. ) 
G(22)=1.5*ALP/(4.*ALP+5.*AL+1. )*(4.*AL+4.*BE+1. ) 
El=l./(SOS+l. ) 

E2=S0S*( (SOS)**3-l. )/( (SOS)**2-l. ) 

E3=Z0/( ( SOS)**2-l. ) 


o n 
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G1=.5*Q41 

G2=.5*Q42 

G3=.5*(Q43+Q151) 

Dl=( (Q41-AN*U41 )-2.*A*(Q31-AN*U31) )/(4.»AL+l . ) 

D6=( (Q42-AN*U42)-2.*A*(Q32-AN»U32) )/(4.*AL+l.) 

D2=( (Q43-AN»U43)-2.*A*{Q33-AN*U33) ) / (4 . *AL+ I . ) + ( (Q151- 

AN»U151)-4.* 

1A*(Q141-AN*U141) )/( i6.*AL+l. » 

D7=-D1*F02  -06*F01**2  -D2*F02*»2 

THE  FOLLOWING  SUBSCRIPTED  VARIABLES  C(IJ,0(I)  ARE  THE 
COEFFICIENTS  OF  EQUATIONS  3.5A  AND  3.5B,  RESPECTIVELY 
C( I ) = .5*( (SOS) **2+1. )*AL/(AL+4. ) 

C(2)=E2*(2.*G( 18)*E3+G( l)*G( 17) )*G(19)/3. 

C( 3)=-2.*E3*G( 14)/CAY**2  +G( 1 ) *Ul 1-G ( 2 ) *U2 I- ( 2 . *E3*G ( 1 

5) /CAY**2  +G 

1( l)*R( l)-G(2)*R(6)+G(3)*R( 1 1 ) +G( 4 ) *R ( 16 ) +G ( 16 ) *G1 ) *F02 

-(G(11*R(2)- 

2G(2)*R(7)+G(3)*R(12)+G(4)*R( 17)+G(16)*G2)*F0l**2-(G(  1) 

*R(3)-G(2)*R 

3(8)+G(3)*R(13)+G(4)*R(18)+G(5)*R(2l)+G(6)*R(23)+G(16)* 

G3)*F02**2 

C(4)=G( 1 )*R(2)-G(2)*R(7)+G(3)*R( 12)+G(4)*R( 17)+G(16)*G 

2 

C( 5 ) = 2.*E3*G( 15)/CAY**2  +G ( 1 ) * ( R ( 1 )-R ( 4 ) +U12 )-G ( 2 ) * ( R ( 

6) -R(9)+U22) 

l+G(3)*(R(ll)+R(14)+U81)+G(4)*(R(16)+R( 19)+U9i ) +G( 16) *G 

1 

C(6)=G(1)*(R{3)+R(5))-G(2)*(R(8)+R(10))+G(3)*(R(13)+R( 

15) )+G(4)»(R 

1( 18)+R(20) )+G(5)*(R(21 )+R(22) )+G(6)*(R(23)+R(24) )+G( 16 

)*G3 

C(7  )=(G( 1)*R(4)-G(2)*R(9)-G(3)*R( 14)-G(4)*R( 19) )*F01-( 

G(1)*R(5)-G( 

12)*R( 10)+G(3)*R( 15)+G(4)*R(20)+G(5)*R(22)+G(6)*R(24) )* 

F01*F02 

C(8)={2.»E3*G( 14)/CAY**2-G( 1 ) »U1 1+G ( 2 ) *U2 1 ) *FO 1- (G ( 1 ) » 

U12-G12)*U22 

1+G( 3)*U81+G(4)*U91 )*F01*F02 
C(9)=IG(1)*G(17)-2.*G( 18)*£3*S0S)*2.*E1/CAY 
D( 1 )=-2./3.*E3*G( 19)*E2-1./3.*G( 1 ) *G { 19 ) *E2  + 9. /32. *G ( 1 

)*AL/(AL+4.) 

1*( ( SOS)**2+l. ) 

D(2)=  (8./12.*G(21)/(4.*AL+9. ) *E3+1 . /6. *G ( 22 ) *G( 19 ) ) *E 

2 

D(3)=-G(20)*£3/(CAY»*2)-G(8)*U31-G(9)*U41+.5*G(8)*U61+ 

.5*G(9)*U71- 

12.*G{7)»G1  +{-G(21 )*E3/(CAY**2)+G( 7)*R(25)+G(8)*R{ 30)+ 

G(9)*R(35)~. 

25*G( 7) *R(70)-.5*G(8)*R(40)-.5*G19)*R{45)+G( 10) *R(50)-G 

( 11)*R(54)-. 

35*G( 10)*R( 58)+.5*G( 1 1 ) *R ( 62 ) -G( 22 ) *G 1 ) *F02  +(-G(7)*R(2 

6)+G(8)*R(31 

4)+G(9)*R( 36)-.5*G( 7)*R(71 )-.5*G(8)*R(41)-.5*G(9)*R(46) 

+G(10)*R(51) 
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5-G(  in*R{55)-.5*G{  10 ) *R(  59 )- . 5*G { 1 1 ) *R  ( 63 ) -G ( 22 ) *G2  ) *F 

01**2  +(-G(7 

6)*R(27)+G(8)*R(32)+G(9)*R(37)+.5*G(7)*R{72)-.5*G(8)*R( 

42)-.5*G(9)* 

7R(47)-G(  10)*R(52)+G(  1 1 ) *R  { 56 ) +.5*G  ( 1 0 ) *R  { 60  )- . 5*G(  ID* 

R(64)+G( 12)* 

8R(66)-G( 13)*R(67)-.5*G( 12)*R(68)+.5*G( 13 ) *R ( 69 )-G( 22 ) * 

G3)*F02**2 

0{4)=G(21)*E3/(CAY*»2)  -G ( 7 ) *R ( 2 5 )-G ( 8 ) * ( R { 30 ) +U33 ) -G ( 

9) *{R( 35)+U4 

13)+.5*G( 7)*(R(70)-U5)+.5*G(8)»(R(40)+U63)+.5*G(9)*{R(4 

5)+U73)-G(10 

2)*(R(50)-U141)*G( 11)*(R(54)-U151 )+.5*G( 10)*(R(58)-U171 

)-.5*&(ll)*{ 

3R(62)-U181)+G(22)*G1-2.*G(7)*G3 
D{ 5)=  G( 7)*R(27)-G(8)*R( 32)-G(9)*R( 37)-.5*G( 7)*R(72)+. 

5*G(8)*R(42) 

1 + .5*G(9)*R(47)+G( 10) *R(52)-G( 1 1 ) *R I 56 ) - . 5*G ( 10 ) *R ( 60 ) ♦ 

.5*G(11)*R(6 

24)-G( 12)*R{66)+G( 13)*R(67)+.5*G( 12 ) *R ( 68 ) - . 5*G ( 13 ) * R ( 6 

9)+G(22) *G3 

0(6  )=  -G(7)*(R(28)-R(26) )-Gl8)*(R(33)+Rl31) )-G(9)*(Rt3 

8)+R(36) )-.5 

1*G( 7)*{R(73)-R(71) ) + . 5*G { 8 ) * ( Rt 43 ) +R (41 ) ) +. 5*G ( 9 ) * ( R ( 4 

8)+R(46) )-G( 

210)*(R(53)+R(51) )-G(ll)*(R(57)-R(55) )-.5*G( 10)*(R(61)- 

R(59) )+.5*G( 

311)*(R(65)+R(63) )+G(22)*G2 

0(7)=(G(7)*R{29)-G(8)*R(34)+G(9)*R(39)+.5»Gl 7)*R(74)-. 

5*G(8)*R(44) 

1+.5»G(9)*R(49) )*F01  + ( G ( 7 ) *R ( 28 ) +G ( 8 ) *R ( 33 ) +G ( 9 ) *R ( 38 ) 

+.5*G(7)*R{ / 

2 3 ) - . 5*G ( 8 ) *R ( 43 ) - . 5*G ( 9 ) *R ( 48 ) *G ( 10 ) *R ( 53 ) +G ( 1 1 ) *R ( 57 ) 

+.5»G( 10)*R( 

361)-.5*G(11)*R(65) )*F01*F02 
0(8)  =-G ( 7 ) *R { 29 ) +G ( 8 ) * ( R ( 34 ) -U32 ) -G ( 9 ) * ( R ( 39 ) +U42 ) - . 5* 

G(7)*(R(74)+ 

1U5)+.5*G(8)*(R(44)+U62)-.5*G(9)* (R(49)-U72)-2.*G(7)*62 
D(9)=  (G(20)*E3/ (CAY**2)  +G ( 8 )*U31+G ( 9 ) *U41-. 5*G ( 8 ) *U6 

1-.5*G(9)*U7 

11+2.*G(7)*G1 )*F02  +(G(8) *U32+G(9)*U42+.5*G( 7)*U5-,5*G( 

8)*U62-.5*G( 

29)*U72+2.*G( 7)*G2)*F01**2  + ( G ( 8 ) *U33+G ( 9) *U43+ . 5*G ( 7 ) * 

U5-.5*G(8) »U 

363“.5*G( 9) *U73-G( 10) *U141+G ( 1 1 ) *U151+. 5»G ( 10 ) *U171-. 5* 

G( 11 )*U181+2 

4.*G(7)*G3)*F02**2 
D( 10)=2.*(2.*E3*SOS-G(7) )*E1/CAY 
D( 1 l)=(G(22)-2.*G(2l )*E3*S0S/AL)*E1/CAY 
THE  SUBSCRIPTED  VARIABLES  Z(I)  ARE  RELATED  TO  THE 
COEFFICIENTS  OF  THE  UNIT  END  SHORTENING  EQUATION  AS 
SHOWN  IN  SUBROUTINE  SOLVE 

Z( 1 )=-(Dl-2.*(l.-AN)*Gl  + (l.  + AN)  /SOS  *AL/CAY**2) 
Z(2)=(1.5*ALP/( 16.*AL+1. )-D2+2.*(l.-AN)*G3) 


o n 


68 


Z( 3 ) = ( .25*AL/(4.*AL  + 1. )* (2. »AL  + l . )-D6+2.*( 1 .-AN)*G2  J 
Z I <»)=-- 1 ( .25*AL/(  4.»AL+l.  )*(2.*AL  + 1.  ) +2 . * ( 1 .-AN )*G2  ) *F0 

1**2  +U.!>*A 

iLP/ (16.*AL+1. )+2.*( l.-AN)*G3)*F02**2+(2.*( l.-AN) *Gl-( 1 

.+AN)/SOS  * 

2AL/CAY**2)*F02  +07) 

Z(5  ) = l./3.*( ( l.-AN)- ( l.+AN) /SOS) *G( IS) 

RETURN 

END 


CSOLVE 

SUBROUTINE  SOLVE 
DIMENSION  CIS), Dill), Z(5) 

COMMON  C,D,Z, T,SOS,RBYH,UO,AN,EM,EN, V,DTAU,F01 ,F02,CAY 

,E1,E2,E3,Z0 

1,CTN 

C PROGRAM  FOR  SOLVING  2 NONLINEAR  SECOND-ORDER  DIFF.  EOS 

C DEFINE  THE  4 FUNCTIONS 

DIFF3FI F1,F2, TAU)=DTAU*1./C1*(C3*FI+C4»F1*«3+C5*F1*F2+ 

C6*Fl*F2**2+ 

iC7*F2+C8+C9*V*TAU*Fl-C2»Fl/(CCl+CC2*F2)*(CC3*F2+CC4*F2 

**2+CC5»F2** 

23+CC6*Fl**2*F2+CC7*Fl+CC8*Fl**2+CC9+{C10+Cll*F2)*V*TAU 

) ) 

DIFF3F(F1,F2,TAU)=DIFFERENTIAL  OF  F3  WHICH  IS  THE 
R.H.S.  OF  EQUATION  (3.8)  (FOR  1=3)  TIMES  DTAU 
DIFF4F(FI,F2,TAU)=DTAU*1./(CC1+CC2*F2)*(CC3*F2+CC4*F2* 

*2+CC5*F2**3 

i+CC6*Fl»*2*F2+CC7*Fl+CC8*Fl»*2+CC9+ (C10+Cll*F2)*V*TAU) 
DIFFIFI VEFl)=DTAU*VEFi 
DIFF2FI VEF2)=DTAU*VEF2 
C1=C( 1) 

C2=C(2) 

C3=C(3) 

C4=C(4) 

C5=C(5) 

C6=C(6) 

C7=C(7) 

C8=CI8) 

C9=C(9) 

CCl=D( 1) 

CC2=D(2) 

CC3=D(3) 

CC4=D(4) 

CC5=D(5) 

CC6=D(6) 

CC7=D(7) 

CC8=D(8) 

CC9=D(9) 

C10=D( 10) 

C11=D( 11) 

EBAR  =2.*CTN**2/CAY 


EBAR=QUANTITY 

VEF1=0 


USED  ro  NORMALIZE  UNIT  END  SHORTENING 


VEFl  = 

VEF2=0 

F1=F01 

F2=F02 

TAU=0 

ZETO=0 


VELOCITY  OF  FI 


I = l 

WRITE  OUTPUT  TAPE  6, 13,T , SOS, RBYH, UO , AN, EM , EN, V, OTAU ,C 

AY,F01,F02,C 

i,D,Z 

13  FORMAT  ( IH1,///(7E16.4//) ) 

WRITE  OUTPUT  TAPE  6,1A 

14  FORMAT  ( IH0,8X,3HTAU,14X,2HF1,14X,2HF2,12X,4HWBYH,12X, 

4HUESN,14X,1 

IHQ//) 

APPLY  RUNGE  KUTTA  TO  SOLVE 
300  P1=DIFFLF( VEFl) 

Q1=DIFF2F( VEF2) 

T1=DIFF3F(F1,F2,TAU) 

U1=DIFF4F(F1,F2, TAU) 

P2=DIFFlF(VEFl+Tl/2. ) 

Q2=DIFF2F( VEF2+U1/2. ) 

T2=DIFF3F(Fl+Pl/2. ,F2*Ql/2. , TAU+DTAU/2. ) 
U2=DIFF4F(Fl+Pl/2.,F2+01/2. ,TAU+DTAU/2. ) 

P3*DIFF1F( VEF1+T2/2. ) 

G3=0IFF2F( VEF2+U2/2. ) 

r3=0IFF3F(Fl+P2/2. ,F2+Q2/2. , TAU+DTAU/2. ) 
U3=DIFF4F(Fl+P2/2.,F2+Q2/2. , TAU+DTAU/2. ) 

P4=DIFF1F( VEF1+T3) 

Q4=DIFF2F(VEF2+U3) 

r4=DIFF3F(Fl+P3,F2+Q3,TAU+0TAU) 

U4=DIFF4F(F1+P3,F2+Q3,TAU+DTAU) 

Fl=Fl+(Pl+2.*P2+2.*P3+P4)/6. 

F2=F2+(Ql+2.*02+2.*Q3+Q4)/6. 

VEFl=VEFl+(Tl+2.*T2+2.*T3+T4)/6. 

VEF2=VEF2+(U1+2.*U2+2.*U3+U4) /6. 

TAU=TAU+DTAU 

UES  = (CTN**2)*(Z( l)*F2+Z(2)*F2**2+Z(3)*Fl*»2+Z(4)+4.» 

E1/CAY*V*TAU 

1-Z( 5)*E2*DIFF4F(F1,F2,TAU)/DTAU) 

UES=UNIT  END  SHORTENING 

Q=V*TAU 

Q = Q-TILDA 

UESN=UES/EBAR 

UESN  = EPSILON  TILDA 

ZET=UESN 

WBYH=(F1+F2)*R8YH*CTN**2 

WRITE  OUTPUT  TAPE  6, 70 , T AU , F 1 , F2 , WB YH, U£ SN , U 
70  FORMAT  (6E16.4) 

IF  (I)  65,50,50 
50  IF  (ZET-ZETO)  55,60,60 
55  I=-l 
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60  ZETO=ZET 

IF(l.5-Q)  1,1,300 
C EXIF  CONTROL 

65  IF  (ZET-ZETO)  60,60,1 

C EXIT  CONTROL 

1 RETURN 

END 


* DATA 

C TYPICAL  DATA  CARDS 

+0.050E+00+1.171E+00+1.750E+02+1.000E-04+0.300E+00+0.800E+01 


+0.100E-01+1.000E+00 


+1.200E+01 
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